Period-doubled breathing in trapped Bose-Einstein condensates 
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The response of a trapped Bose-Einstein condensed gas to a periodic driving force is studied 
theoretically in the framework of the nonlinear Gross-Pitaevskii equation. The monopole mode 
is driven by periodical modulation of the frequency of the isotropic harmonic trapping potential. 
Period-doubled motion is shown to occur when the driving is strong and close to resonance with the 
monopole mode. A variational model predicts chaotic oscillations but is found to disagree with the 
full solutions to the Gross-Pitaevskii equation. 
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A trapped Bose-Einstein condensed gas is a damped, 
nonlinear oscillator. The modes of oscillation result 
from an interplay between dispersion, the repulsive inter- 
particle interactions and the confining trap potential Q, 
and the damping is offered by friction against thermal 
atoms 0, Q • Interesting physics can therefore be antici- 
pated if the trap is modulated so as to drive the modes 
of oscillation: periodically forced, nonlinear oscillators 
can generally be expected to exhibit period-doubled mo- 
tion and chaos for certain parameter values. Textbook 
examples include the Duffing oscillator and the forced, 
damped pendulum £|. 

In this paper, we shall study the oscillations of Bose- 
Einstein condensed gases subject to periodic driving 
forces, with particular attention to the conditions for 
period-doubled motion and the onset of chaos. The ba- 
sis for all calculations is the Gross-Pitaevskii equation 
(GPE) for the condensate wave function ip(r,t), 
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The wave function ip(r, t) determines the local fluid den- 
sity and velocity through n(r, t) = \ip(r, t)\ 2 and v(r, t) = 
(7i/m)V arg-0(r, t); it is normalized to the number of con- 
densed atoms N. The nonlinear term in the GPE de- 
scribes inter-particle interactions in the s-wave approxi- 
mation and its strength is determined by the s-wave scat- 
tering length a. The trap potential V(r, t) is a harmonic- 
oscillator potential which in this study shall be taken to 
be isotropic and periodically modulated in time: 



V(r, t) = ^muj 2 r 2 (1 + A cos tit) , 
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Assuming that the density profile stays fixed and only its 
overall radial extent is allowed to vary, the cloud size is 
described by an effective equation Q, 



R= -R(l + A cos fit) 



R 4 



2-fR. 



(3) 



This equation is often derived in a variational scheme Q 
and will therefore be referred to as "variational" in the 



following. The first term on the right hand side repre- 
sents the trap potential including the periodic modula- 
tion. The inverse-power term originates from the inter- 
particle interactions. The kinetic energy associated with 
the density gradient is assumed negligible in compari- 
son with the trap and interaction energies and the cor- 
responding term is therefore not included in Eq. 
this approximation holds when the scattering length is 
large and the mean density is high |§j]. The unit of 
length is taken to be the Thomas-Fermi radius Rtf — 
a osc (15iVa/a osc ) 1 / 5 , where a osc — (h/muj) 1 ^ 2 is the oscil- 
lator length, and the time and frequency have been scaled 
by the trap frequency u>, so that ti — ti/u> and t = cut. 
We shall hereafter work only in these dimensionless units, 
and we therefore immediately drop the bar on t and ti. 
The final linear friction term in Eq. is introduced by 
hand and represents the dissipation due to the thermal 
gas. For the physical situation at hand, damping plays a 
pivotal role and it is therefore important to allow for this 
effect. This phenomenological treatment of the dissipa- 
tion is known to well reproduce the experimental findings 
of Ref . |2( . The friction is provided by Landau damping 
and depends on temperature through the relation 
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where uj m is the frequency of the mode in question, and T c 
is the critical Bose-Einstein condensation temperature. 
The damping coefficient 7 = 0.02, which is used in much 
of the calculations of this paper, corresponds to a cloud 
of N — 10 6 atoms with a osc /a — 100 at a temperature 
T = 0.1T C . 

The virtue of the variational equation © is of course 
its simplicity, allowing for fast numerical solution, so that 
large parameter spaces can be covered, and in a few in- 
stances the existence of analytical solutions. However, 
the results obtained this way must of course be double- 
checked against the solution of the full Gross-Pitaevskii 
equation, Eq. l(T]l. to find out whether deformation of 
the density profile plays a crucial role. For the present 
case of isotropic forcing, we shall assume that the cloud 
stays spherically symmetric. The Gross-Pitaevskii equa- 
tion, thus made effectively one-dimensional, is propa- 
gated in time with the Crank-Nicolson method. Dissi- 
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pation is represented by letting the time be complex, so 
that in solving the discretized equation, the time step is 
At = A<tj(1 + iTgp)- It is found numerically that the 
parameter Tgp is related to the damping of the breathing 
mode through 7 sw 3Fgp- 

The problem depends on three physical parameters: 
the amplitude A of the perturbation, the frequency f2 of 
the perturbation, and the damping 7 0. Luckily, the 
dynamics of the cloud is found not to depend sensitively 
on the latter, which is difficult to control or even mea- 
sure with great precision. When the system is not driven, 
A = 0, the system performs damped sinusoidal oscilla- 
tions with the frequency of the breathing mode, which is 
equal to y/E in the strong-coupling limit considered here. 
When the amplitude of the driving force A is nonzero, 
but not too large, we expect the motion to be frequency 
locked, i. e. it performs breathing motion that is peri- 
odic with the frequency O, or possibly fractions thereof, 
in which case we have period-multiplied motion. When 
A > 1, the instantaneous trapping potential will for a 
certain time interval during each period not be confin- 
ing, because the instantaneous squared trap frequency, 
ui 2 (l + Acosfii), turns negative. This does not, how- 
ever, necessarily imply that atoms are lost, because such 
events happen on a finite time scale: if the trap is opened 
but then sufficiently rapidly closed again, the atoms re- 
main confined. For large enough A, the system does 
indeed turn unstable, but the threshold value depends 
sensitively on f2 (and less sensitively on 7). 

The region of instabilityof the variational Eq. © have 
been examined in Refs. |ToL ITU]. For small forcing am- 
plitudes A, the instabilities occur in the vicinity of the 
frequency values f2 = 2/n, n = 1,2,3..., and as A is 
increased the instability regions grow larger. Since insta- 
bility is associated with large-amplitude motion, the rel- 
evant frequency is in fact the bare trap frequency (which 
is 1 in the units currently employed), rather than the 
quadrupole frequency y/E, and the instability regime is 
the same as that for the Mathieu equation. The inclusion 
of damping somewhat shrinks the instability regimes. 
This analysis predicts an appreciable domain for stability 
even for values of A well above unity, especially when Q is 
large. In fact, when f2 approaches infinity, the threshold 
value of A for instability also approaches infinity. 

Figure [T] shows examples of the time evolution of R(t), 
the solution of Eq. ©, for two instances of the values 
of the parameters O and A. The limit cycle, i. e. the 
steady orbit that is approached at long times, in the 
phase space spanned by R and R is also shown. The 
damping parameter is fixed at 7 = 0.02. In the first ex- 
ample, the forcing frequency = 2.2 is comparable to 
the quadrupole frequency and the forcing amplitude is 
rather low, and therefore the motion is in this case nearly 
sinusoidal (so that the limit cycle takes on a nearly el- 
liptic shape) and locked to the forcing frequency. For a 
larger amplitude A the shape of the oscillations is dis- 
torted due to the strongly anharmonic effective potential 
well that R resides in. When the forcing frequency is 
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FIG. 1: Time evolution and limit cycles of the forced conden- 
sate, with driving amplitude A and frequency Q as indicated 
above the panels. The leftmost plots display a time span of 
eight periods of the driving force. 



smaller than an integer fraction of the breathing-mode 
frequency, SI < y/E/n, an n-th harmonic is seen if A is 
large enough (but the motion is still periodic with period 
fi). This is in Fig. □ illustrated for (1 = 0.53 = \/5/4.2 
where four loops are seen in the limit cycle. These varia- 
tional results have been compared to the solution of the 
full GPE 0}, and are in excellent agreement. 

For larger values of the forcing amplitude A more com- 
plex behavior is observed. When f2 is close to reso- 
nance with the breathing-mode frequency ft = y/E or 
fractions thereof 01 ■> an d A is large enough, the varia- 
tional system @ exercises period-doubling bifurcations 
and chaotic motion. Figure [5] shows the period-2 mo- 
tion seen at Q = 1.25 w VE/2, 7 = 0.02 and A = 1.5. 
Shown are both the variational result for R and the rms 
radius, yj (r 2 ), of the cloud as calculated using the Gross- 
Pitaevskii equation, scaled by the Thomas-Fermi value 
Rrms = -\/3/7i?TF to make the scale agree with that of 
the variational study. The limit cycle is plotted in the 
space of R(t) and R(t — At), where At is a small time 
delay; this yields an image tilted by 45 degrees compared 
to that in R-R space. The variational model still agrees 
well with the full GPE. However, concerning chaos the 
situation is different: the solutions to the full GPE fail to 
exhibit chaotic motion when the variational does. An ex- 
ample is shown in Fig. where the parameters f2 = 0.7, 
A = 1.8 and 7 = 0.02 are chosen. That the solution to 
the effective-radius equation is indeed chaotic has been 
checked by calculating the Lyapunov exponent. How- 
ever, the solution to the full Gross-Pitaevskii equation is 
seen to be perfectly regular: still period-doubled, albeit 
with a rather large amplitude. The grid parameters have 
been chosen so that the wave function both in configu- 
ration and momentum space stays well inside the grid 
boundaries at all times; grids up to 1500 points with a 
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FIG. 2: Time evolution and limit cycles of the forced conden- 
sate. Variational results are shown in the upper panels, while 
the lower panels show the results of the full Gross-Pitaevskii 
equation. The delay parameter for the limit cycle in the lower 
plot is chosen as At = 0.5. 
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FIG. 3: Time evolution and limit cycles of the forced conden- 
sate. Variational results are shown in the upper panels, while 
the lower panels show the results of the full Gross-Pitaevskii 
equation. 



grid constant Ax = 0.04a osc have been tried. 

The result is quite general: nowhere in phase space 
does chaotic motion, or even periodic motion with period 
longer than 2, occur in the GPE solution. The discrep- 
ancy is not due to the omission of the dispersion-related 
kinetic-energy term in Eq. (J3J; that is readily checked by 
including such a term whereby the solution still exhibits 
chaotic motion. Rather, the fallacy of the variational Eq. 
© is that it does not let the wave function change its 
shape during the dramatic changes of the trap potential 
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FIG. 4: Density distribution at a few time instances, as calcu- 
lated by the full Gross-Pitaevskii equation. Radial real-space 
density distributions are shown to the left and momentum- 
space distributions to the right. Topmost panels: periodic 
motion with — 0.71 and A = 1.3. Lower panels: period-two 
motion with Q — 0.71 and A = 1.6; this case was predicted by 
the variational equation to be chaotic. The different curves in 
each panel are chosen at arbitrary equally spaced instances, 
but the same set of time instances were used in all panels. 

between confining and repelling. This may come as a sur- 
prise, because a variational study is almost always seen 
to faithfully mimic the physics at least qualitatively 0. 
A closer look at the wave function gives an explanation. 
Figure 0] shows a few snapshots of the wave function for 
the cases of period-one and period-two motion. In the 
regime where there is agreement between the variational 
equation and the GPE, the profile of the wave function 
does not change drastically, in accordance with the as- 
sumptions behind the variational equation, but in the 
regime predicted to be chaotic, waves propagate back 
and forth through the cloud. Clearly the variational sys- 
tem has too few degrees of freedom to faithfully model 
the behavior of the cloud in this regime. In a manner 
of speaking, the instability that manifests itself in short- 
wavelength density waves in the full GPE study, must 
instead find its outlet as chaotic motion in the restricted 
variational study. 

We conclude by showing in Fig. the full phase dia- 
gram of the forced system as calculated by the variational 
equation. The damping parameter is still 7 = 0.02. Small 
regions of period-multiplied motion and chaos are seen to 
exist around the stability tongues just above f2 = yE/n 
for all integers n > 2. The largest chaotic regime is also 
indicated although it was found to be an artifact of the 
variational equation. 

To summarize, the Gross-Pitaevskii equation in a trap 
with a periodic modulation of the trapping potential has 
been seen to support period-doubled oscillations, but 
chaotic motion was not observed. In contrast, a sim- 
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FIG. 5: Phase diagram for the forced condensate. Solid lines 
enclose the regime of stable motion; dashed lines enclose the 
regimes of period-doubled motion and the dotted line encloses 
the regime where the variational study displays chaotic mo- 
tion. The insets labeled (a) and (b) zoom in on two important 
regions. 

plified description of the same system using a collective 
coordinate predicts chaotic motion and period- multiplied 
motion with period higher than two, because of its failure 
to describe short- wavelength density modulations. In- 
deed, it was suggested in Refs. 0, that chaotic mo- 



tion cannot occur in Bose-Einstein condensed systems; 
however, it seems that the reasoning leading to that con- 
clusion is not applicable in the present case. Reference 
|13| proposed that the exponentially growing separation 
between neighbouring orbits in phase space that is char- 
acteristic of chaos implies that the depletion of the con- 
densate grows equally fast, and Ref. 11] suggested that 
this can be taken as proof that chaotic oscillation in Bose- 
Einstein condensates is impossible. However, said anal- 
ysis does not address the case of bounded chaotic oscil- 
lations where the exponential separation of nearby or- 
bits only takes place over short times while at long times 
the whole bundle of orbits stays bounded, as in Fig. 
Furthermore, if the excited modes are coherent with the 
condensate, it is not at all clear that their growth must 
be identified with depletion of the condensate. The ab- 
sence of chaos observed here should therefore be seen as 
an accidental fact, and the question of whether chaotic 
oscillations in fact can be observed for slightly different 
physical situations - for example exploring other param- 
eter values or driving modes with higher multipolarity - 
still remains open. 
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